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0.  Introduction 

— ^The  traditional  task  of  geodetic  astronomy  as  seen  from  the  viewpoint  of  prac¬ 
tical  geodesy  is  the  determination  of  three  spatial  orientation  parameters  of 
a  vertically  set  up  observation  instrument  relative  to  a  global  reference  frame 
fixed  in  the  earth.  For  this  horizontal  and  vertical  directions  are  measured 
at  registered  instants  tc  stars*  whose  coordinates  are  assumed  to  be  known  in 
a  space  fixed  reference  frame,  the  time  dependent  orientation  of  which,  rela¬ 
tive  to  the  earth  fixed  frame,  is  also  presumed  to  be  known. 


~*The  required  parameters  are  the  Jttnonomic al  longitude.  A' and  the  ou6XA.unoiM.caZ 
latitude-#*  which  fix  the  direction  of  the  local  gravity  vector  relative  to  the 
earth  fixed  reference  frame,  and  the  horizontal  orientation  unknown  £  of  the  in¬ 
strument,  which  yields  the  astronomical  azimuths  A  in  connection  with  the  mea¬ 
sured  horizontal  directions  to  terrestrial  objects. 

A  spatial  reference  frame  is  here  understood  as  a  triad  of  orthonormal  base  vec¬ 
tors,  which  is  fixed  to  a  distinct  origin  point  and  which  is  taken  as  being  ro- 
tationally  and  translationally  invariable  in  time. 

The  models  used  traditionally  in  geodetic  astronomy  are  described  in  the  stan¬ 
dard  textbooks  (for  example,  K.  Ramsayer,  1970,  see  Appendix  A.1)  mainly  by 
making  use  of  spherical  trigonometry.  >The  presentation  used  in  this  report  is 
based  upon  orthonormal  triads  of  base  vectors  of  different  reference  frames.  A 
compact  and  strictly  systematical  presentation  is  obtained  by  a  commutative  dia¬ 
gram  of  transformations  between  the  respective  bases>  which  is  introduced  in 
chapter  1  .^Fundamental  relations  between  the  parameters  of  geodetic  astronomy 
result  from  this  commutative  diagram.  "In  chapter  2  these  relations  are  trans¬ 
formed  iff ter  the  linearization  with  suitable  approximate  values ^fn  such  a  way 
that  a  system  of  condition  equations  with  unknowns  ensures,  the  Gau&-HelmeAt 
model  of  the  adjustment  calculation.  ^This  system  is  specified  for  different 
combinations  of  observations.  In  chapter  3  the  accuracies  in  the  determination 
of  A,4>  and  Z,  which  are  expected  for  different  configurations  of  stars,  are  es¬ 
timated  in  simulation  studies  according  to  the  GauB-Helmert  model  and  then  the 
results  are  illustrated  in  diagrams. 
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The  underlying  principles  for  this  report  were  taken  from  the  dissertation  of 
8.  ZichteA,  "Entwurf  eines  nichtrelativistischen  geodatisch-astronomischen 
Bezugssystems",  Deutsche  Geodatische  Kommission  Heft  C322,  Munich,  1986,  and 
from  the  manuscript  of  the  lecture  "Geodetic  Astronomy",  which  B.  Richter 
gives  at  the  University  of  Stuttgart,  Federal  Republic  of  Germany. 
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1 .  Relations  between  the  reference  frames  in  geodetic  astronomy 

First  of  all  in  this  chapter  the  fundamental  relations  in  geodetic  astronomy 
between  observations,  unknowns  and  given  coordinates  shall  be  derived  which 
will  be  needed  further. 

1.1  The  systematical  structure  of  the  reference  frames 


The  reference  frames  used  in  geodetic  astronomy  may  be  arranged  on  different 
levels  which  are  numbered  in  turn  or  indicated  by  symbols:  0  corresponds  to  ' ,  Xj 

1  to  *,  2  to  *,  3  to  o.  One  fundamental  vector  V1  belongs  to  every  level  i.  n 

In  details  this  is  as  follows: 
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o 


V' 


=  Z 


V 


1 


v*  = 


r 


V 


2 


=  V* 


= 


V 
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the  position  vector  from  the  point  of  observation  to 
the  target  object  (terrestrial  or  celestial),  which  is 
generally  a  star; 

the  negative  gravity  vector; 

the  earth  rotation  vector  (it  has  the  direction  of  the 
axis  of  the  earth,  points  to  the  North  Pole  and  has  the 
value  of  the  earth  rotation  rate); 

the  ecliptic  normal  vector  (it  points  to  the  northern 
pole  of  the  eel iptic). 


An  orthonormal  reference  frame  E1  belongs  to  every  level  with  its  base  vectors 
as  follows: 


~3 


=  norm  V1 


S 


norm  (V1+1  x  V1) 


ti 


E1  x  E1 
~2  *  t3 


(1-D 

(1-2) 

(1-3) 


Here  "norm"  denotes  the  abbreviation  for  the  normalization  of  a  vector,  and  "x" 
the  vector  product. 


New  reference  frames  are  at  the  lower  end  the  observational  frame  E'  of  the  level 


* 


"0",  whose  third  base  vector  is  located  in  the  direction  of  the  observation  and 
which  is  unique  since  it  is  not  a  reference  frame  in  the  literal  sense,  because 
there  are  no  vectors  described  with  regard  to  this  frame,  and  at  the  upper  end 
the  ecliptic  frame  E°,  which  has  hardly  any  practical  importance. 

In  addition  to  the  systematical  E-triads,  F-triads  also  appear  on  each  level. 
These  systems  have  the  common  third  base  vector  with  the  appertaining  E-frame, 
nevertheless  the  direction  of  the  first  and  second  base  vector  does  not  follow 
from  the  systematic  structure  of  the  fundamental  vectors  V1,  V1+1 ,  but  from  a 
more  or  less  arbitrary  definition. 


A  new  F-triad  is  the  theodolite  frame  F*,  whose  first  base  vector  F, .  lies  in 
the  direction  "zero"  of  the  azimuth  circle  of  a  theodolite  which  is  set  up  in 
the  astronomical  horizon.  The  longitudinal  angle  (see  below)  of  an  observed 
direction  in  the  local  horizon  frame  is  the  horizontal  direction  T  and  is  re¬ 
corded  systematically  -in  the  clockwise  direction,  but  conventionally  courvteA- 
clockwise,  T  =  -T  .  The  latitude  angle  (see  below)  is  the  vertical  direction 

S  C 

as  in  the  horizon  frame  E*. 

AC 


The  transformation  from  a  frame  E1  to  the  appertaining  frame  F1  is  always  a 

at  a 

counter-clockwise  rotation  round  the  common  third  axis  with  the  orientation 
angle  H1, 


F1  =  R3(Hi)Ei  . 


(1-4) 


R3  is  the  rotation  matrix,  which  describes  a  rotation  of  a  frame  round  the 
third  axis.  It  is 


R3(y) 


cosy  siny  0 

-siny  cosy  o 

0  0  1 


(1-5) 


Corresponding  to  eqn.  (1-5)  the  rotation  matrices  for  the  rotations  round  the 
first  and  second  axis  read 


■yMBGMOOf*: 


,  R2  and  g3  are  also  called  elementary  rotations.  The  orientation  angles  H 
(see  eqn.  (1-4)  )  are  in  detail: 


H1  =  H*  =  E  the  orientation  unknown  of  the  theodolite  which  has 

been  set  up; 

H2  =  H  =  the  Greenwich  sidereal  time; 

Gr  ,s 

H3  =  H°  the  angle  between  the  line  of  intersection  of  the  eclip¬ 

tic  with  the  mean  galaxy  plane  and  the  direction  to  the 
vernal  equinox  ±90°. 

For  the  transformation  from  a  frame  E1+1  to  the  underlying  frame  E1  one  needs 

.  —  .  ac 

the  longitudinal  angle  x^+1  and  the  latitude  angle  <^+1  of  the  fundamental 

vector  V1  with  regard  to  the  frame  E1+1: 


I1  -  Rj(90--*itl)  R3(x^)  fi+1 


-i+1 


i+l 


(1-8 


R^,  is  the  special  case  of  a  rotation  matrix  of  Eulerian  type,  in  which  the  three 
elementary  rotations  are  connected  in  a  row  as  follows: 

-  first  rotation  round  the  third  axis,  R3(y  ) 

-  second  rotation  round  the  new  second  axis  with  the  angle  (90°-8)  ,R2(90°-B) 


third  rotation  round  the  new  third  axis,  R  (y  )  • 


The  £  and  F  reference  frames 


I  HFJW.KFl 
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In  the  matrix  R  the  second  rotation  round  the  second  axis  would  take  place 

-e  r 

with  the  angle  B,R_(B). 


These  longitudinal  and  latitude  angles  are  in  detail: 


X?  -  x;  =  As 


-  n  -  b 


$*  =  $ 


the  azimuth  of  the  observational  direction; 


the  vertical  direction 


the  sidereal  time 


the  astronomical  latitude. 


For  the  transformation  from  a  frame  Ei+1  to  the  frame  F1  lying  diagonally  under- 

st  . 

neath,  one  needs  additionally  the  orientation  angle  H1  (see  Fig.  1): 


r  -  -RE(xi+1.^1.Hi)ii*1 


(1-9) 


For  the  transformation  from  a  frame  F1+1  to  the  frame  E1  lying  diagonally  under¬ 
neath,  one  needs  the  longitudinal  angle  A^+1  and  the  latitude  angle  $*+i  of  the 
fundamental  vector  V1  with  regard  to  the  frame  F1+1  : 


I1  - 


(1-10) 


The  latitude  angles  are  the  same  as  above,  the  longitudinal  angles  are  in  de¬ 
tail  : 


A.u  =  a;  =  T 

1  *  s 


the  horizontal  direction  of  the  observation  direction, 
systematically  measured  counter-clockwise,  conventionally 
in  the  clockwise  direction,  T  =  -T  ; 

S  C 

the  astronomical  longitude; 


=  90°-e 


the  orthogonal  complement  to  the  inclination  of  the 
eel iptic. 


E4 


level  4,  mean  galaxy  frame 


level  3  (°),  ecliptical  frame 


level  2  (°),  equatorial  frame  (fixed  in 
space,  fixed  in  the  earth) 


level  1  ('),  horizontal  and  theodolite 

frame 


1  ;  Commutative  diagram  with  reference  frames  in  geodetic  astronomy 


For  the  transformation  from  a  frame  F1+1  to  the  underlying  frame  F1,  one 
needs  additionally  the  orientation  angle  H1: 


c1  • 


(1-id 


The  transformation  in  the  opposite  direction  results  from  the  respective  trans¬ 
posed  of  the  rotation  matrix. 


1.2  The  commutative  diagram  and  the  fundamental  equation  of  geodetic  astronomv 


In  the  rotations  introduced  up  to  this  point  the  right  ascension  a  and  the  de¬ 
clination  6  are  still  missing;  these  describe  the  observation  direction  E3, 
to  a  star  with  regard  to  the  space  fixed  equator  frame,  namely  in  the  same  way 
as  A  and  B  do  this  with  regard  to  the  horizontal  frame.  Indeed,  the  rotation 
R^a.SjE*  does  not  lead  to  the  frame  E';  in  E" =  RE(a,6)E*  the  second  base  vec¬ 
tor  E2„  lies  in  the  equatorial  plane,  but  in  E‘  the  vector  E2,  lies  in  the  hori¬ 
zontal  plane  (east).  Of  course  E'  and  E"  have  the  common  third  base  /ector, 
but  they  differ  by  a  rotation  round  this  vector  with  an  angle  E,  : 


E"  -  R3U)E* 


(1-12) 


A  serial  connection  of  several  transformations  is  called  a  diagram  in  the  alge¬ 
bra.  If  the  transformations  which  have  to  be  reversibly  unequivocal,  form  a 
closed  circle,  this  is  called  a  commutative  diagram.  Then  one  is  able  to  express 
a  transformation  by  means  of  the  others.  Such  a  commutative  diagram  is  presen¬ 
ted  in  Fig.  1  by  the  lines  which  are  thickly  marked.  For  example  the  trans¬ 
formation  E*  -*■  E'  can  be  expressed  by 

E'  =  Re(T,B)R3(Z)E*  (1-13) 

or 

V  =  R3(C)RE(a,6)R3(eGr)RE(A,<i>)E*  .  (1-14) 

As  the  representation  of  the  triad  E'  with  regard  to  the  triad  E*  is  unequi¬ 
vocal,  it  follows  that; 


(1-15) 


RB(TiB)R3(Z)-R3(5)RB(a,6)R3(e  )g(A,*) 


This  equation  reads  with  the  right-hand  side  written  in  full  length 


RE(T,B)R3(l)=R3(-C)R2(90o-6)R3(a)R3(-9Gr)R3(-A)R2(<t-90°) 


(1-16) 


and  can  again  be  reduced  to 


RE(T+Z,B)=R^(C,6-90o)Bg(6G  +A-a,$) 


(1-17) 


These  are  the  desired  fundamental  AeZcutiovu  between  the  parameters  appearing 
in  geodetic  astronomy. 


V//'' 
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2.  The  observation  equations  of  geodetic  astronomy 

The  fundamental  equation  (1-17)  consists  as  a  matrix  equation  of  nine  sepa¬ 
rate  equations,  of  which  only  three  are  independent  of  each  other  because  of 
the  property  of  orthonormality  of  the  rotation  matrices.  These  three  inde¬ 
pendent  equations  represent  condition  equation*  with  unknown*  for  every  star 
(if  T,B  and  8_  are  measured  at  one  instant). 

The  matrices  on  the  left  and  right-hand  side  of  equation  (1-17)  read  as  fol¬ 
lows  when  they  are  multiplied  respectively: 


cos(l+T)sinB 


-sin(I+T) 


cos(I+T)cosB 


sin(E+T)sinB 

cos(E+T) 


sin(I+T)cosB 


-cosB 


Column  1 : 


sin$cos(6Gr+/  :x)sin6cos£-sinii|sin(8Gr+A-a)sin£+cos(J>cos6cos£ 
sin$cos(6„  +A-a)sin6sinC+sin$sin(0^  +A-a)cos£+cos<J>cos<5sin£ 
s  i  n$cos  ( 8Gr+A-a)  cosS-cosOs  i  n<5 

Column  2: 


-sin(0Gr+A-a)sinficos£-cos(8Gr+A-a)sin£ 
-sin(8Gr+A-a)sin«Ssin£+cos(  0Gr+A-a)cos£ 
-sin(0Gr+A-a)cos<5 


Column  3: 


cos$cos(8  +A-a)sin5cos£-cos<J>sin(8„  +A-a)sin£-sin$cos6cos£ 

\j1T  CjJT 

cos<I>cos(8Gr+A-a)sin6sin£+cos<1|sin(0Gr  +A-a)cos£-sin$cosSsin£ 
cos<t>cos(8Gr+A-a)cos<5+sin<i>sin5 


In  the  right-hand  matrix  the  elements  in  the  third  row  are  the  shortest  and  at 
the  same  time  the  only  ones  which  do  not  contain  the  angle  £.  Therefore,  it 
is  the  obvious  thing  to  do  to  select  two  equations  from  this  row  as  indepen¬ 
dent  equations.  As  a  third  equation  one  could  take  an  element  from  another 
row  of  the  matrices  whereby  the  angle  £  ,  in  which  one  is  not  actually  in¬ 
terested,  would  indeed  appear  as  an  additional  unknown.  So  one,  therefore, 
dispenses  with  such  an  equation  and  there  remain  only  two  independent  equa¬ 
tions  for  one  complete  observation  (T,B  and  e_  )  with  the  three  unknowns  A, 

Gr 

sinB  =  cos$cos(  0_  +A-a)cos6+sin<tsin6 

Gr 

sin(I+T)cosB  =  -sin(0^  +A-a)cos6 

Gr 

cos(I+T)cosB  =  sin4>cos(0„  +A-a)cos6-cos<J>sin6 

Gr 


(2-1) 

(2-2) 

(2-3) 


Equations  (2-1)  and  (2-2)  are  independent  of  each  other,  equation  (2-3)  is 
dependent  on  them  both.  It  will  be  used  later  only  for  the  determination  of 
approximate  values.  The  appearing  variables  be  summarized  once  more: 


A 

$ 


astronomical  longitude 
astronomical  latitude 
right  ascension  of  the  star 
declination  of  the  star 
hour  angle 

orientation  unknown  of  the  instrument  (theodolite) 

horizontal  direction;  observed 

azimuth 

vertical  direction,  angle  between  horizon  and  star;  observed 
Greenwich  apparent  sidereal  time;  observed 


2.1  Linearization  and  matrix  representation 

The  equations  (2-1)  and  (2-2)  are  now  linearized  by  developing  them  into  a 
Taylor  progression  (and  stopping  after  the  first  order  term): 


a 


TSWWV’ 
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sinB0+cosBQ6B  4  cos<t>0cosh0cos6+sin$0sin6 

-cos<J>0sinh0cos6(66Gr+6A)  (2-4) 

+( -s  i  n<t>0cosh0cos5+cos4>0  s  i  n<5 )  64> 

sin(S  +T  )cosB.+cosAr.cosBr.(6I+6T)-sinA.sinBri5B 

- - 2.  oo  oo  (2-5) 

4  -sinh„cos6-cosh„cos6(66„  +6A) 

0  0  Gr 


Equations  (2-4)  and  (2-5)  contain  the  terms  of  (2-1)  and  (2-2)  taken  at  the 
point  of  developing  (underlined).  By  assuming  that  there  are  given  approxi¬ 
mate  values  for  the  unknowns  A,0  and  Z  (these  are  A0,<t>0,Z0)  one  can  calculate 
values  for  B  and  T,  so  that  equations (2-1 )  and  (2-2)  (and  (2-3)  because  of  the 
ambiguity  of  sine  and  cosine)  are  satisfied  for  these  approximate  values.  There¬ 
fore,  in  equations  (2-4)  and  (2-5)  the  terms  underlined,  taken  at  the  Taylor- 
point,  compensate.  So  it  follows  from  (2-4)  and  (2-5)  if  one  in  addition  places 
equation  (2-3)  in  (2-4): 

-cosBo6B-cos$0sinh0cos<5(<5eGr+<5A)-cosA0cosB064>  =  0  (2-6) 

-cosA_cosB  (6Z+6T)+sinA„sinB/,6B-cosh„cos5(60„  +6A)  =  0  (2-7) 

u  u  u  u  0  or 

From  these  equations  one  obtains  the  onset  for  the  conditional  equations  in 
respect  of  the  adjustment  problem  by  introducing  the  vector  e  of  inconsistency. 
6B,  6T  and  60Gr  represent  in  this  case  the  difference  between  the  actual  ob¬ 
servations  to  one  star  at  one  instant  and  the  approximate  values  which  go  into 
the  coefficients  of  the  conditional  equations: 


■* 

6B  =  B  -  EL 

0 

constitute  the  vec ton.  o&  obteAvatLoru  y  = 

-6T 

5T  *  T-To 

=>  E{y}  =  y  -  e  . 

-6B 

68=r=  8Gr  "  8GrO 

-66 

Gr 

This  is  taken  into  consideration  in  equations  (2-6)  and  (2-7);  these  equations 
are  now  remodelled: 


(2-8) 


! 


,  Mw'tb  •*_ 
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-cos4>0sinh0cos6»6A-cosA0cosB0*6$-cosB0*eB-cos$0sinh0cos6*eeGr 


=  +cosB  •6B+cos<i’  sinh.cos6*60_ 

0  0  0  Gr 


-cosh0cos6*6A-cosA0cosB0,6Z-cosA0cosB0*eT+sinA0sinB0*eE 


(2-9) 


-cosh.cosS'e^  =  cosA  cosB  •  6T-s i nA  s i nB .  •  <$B+coshrt  cos<5  •  60„ 

0  tfOr  u  u  0  0  0  Gr 


If  observations  are  carried  out  in  respect  of  several  stars  the  general  re¬ 
presentation  for  the  GauR-Helmext  modeZ  of  conditional  equations  with  unknowns 
is 


Ax  +  Be  =  By  -  c  (2-10) 

Equation  (2-10)  runs  in  the  case  of  consistent  approximate  values  (c  =  0): 

Ax  +  Be  =  By  :  B*E{y}  =  Ax  (2-11) 


A  and  B  are  the  coefficient  matrices  for  all  observations  to  all  stars  obser¬ 
ved.  In  the  case  of  three  unknowns  A,4>,I  and  of  observations  to  n  stars  A 
has  the  size  2n  x  3  and  B  2n  x  3n  and  contain  submatrices  Ai  and  B±  for  every 


triplet  of  observations  to  every  star: 


(2-12) 


"B.  ! 

2x3  ! 

1 

1 

"V; 

1 

2x3  : 

!  B3  i 

!  2x3 ! 

i 

1 

i 

(2-13) 


aw 


wwi 


The  elements  a.,  of  the  submatrices  A.  and  the  elements  b„  of  the  submatrices 

ik  3  x-m 

Bj  respectively  are  now  introduced. 

In  schematized  form  it,  therefore,  follows  that: 


au«A«126i.bi2€B*b13e9Gr  =  -b12«B-b136eSr 


a216AM236£+b2ieTtb22£Btb23£0Gr  =  •b216T-b22iB'b2360Gr 


In  matrix  form  this  then  becomes: 


a  a  0 
11  12 


a2i  0  a23 


0  b  b 

12  13 


b21  b22  b23 


0  b12  bl 3 


b21  b22  b23 


-6B  (2- 


No  summation  over  the  same  indices! 


Hence  AJ  and  B.  are 
j  3 


-cos<&„sinh  cos6  -sin<t  cosh  cosiS+cosi>  sin6 
0  0  0  0  0 


-coshQcos6 


-cosAqcosB0 


■cosAocosBq 


-cosB, 


sinAQsinB0 


-cos$0sinh0cos6 


-coshQcos6 


One  perceives  that  an  =  b13,  a  }  =  b23,  a23  =  b 


2.2  The  different  types  of  observations 
One  can  observe  for  eve/cy  star 


a)  vertical  direction  and  time  of  observation  or 

b)  horizontal  direction  and  time  of  observation  or 

c)  horizontal  direction,  vertical  direction  and  time  of  observation  or 

d)  horizontal  direction  and  vertical  direction  (without  time). 

Case  c)  corresponds  to  the  construction  in  (2-14)  to  (2-19).  The  other  cases 
are  computed  as  follows.  Thereby  the  starting  point  was  taken  as  equation  (2-6) 
=  I  and  equation  (2-7)  =  II,  because  one  is  only  interested  in  the  respective 
matrices  A.  and  B.. 


2.2.1  Case  a):  Vertical  direction  and  time 

For  every  pair  of  observables  (that  is  the  vertical  direction  B  and  the  side¬ 
real  time  0Gr)  one  obtains  one  conditional  equation  of  type  I: 

sinB  =  cos$cos(0Gr+A-a)coS(5+sin$sin6  (2-20) 

which  gives  in  linearized  form 

-cosB,6B-cos<t>  sinh  cos<5(60  +SA)+(-sin<i>  cosh  cos5+cos$  sin6)6$  =  0  .  (2-21) 

u  u  o  t>r  OOO 

Then  the  matrices  read 


B.  =  Cb  b  ]  =  C-cosB  -cos<i>  sin(8„  +A  -a)cosSJ 
3  12  13  0  GrQ  0 

Aj  =  Cail  ai2'3 

=  C-cos$  sin(e  +A  -a)cos6  -sin$ncos(e_  +A0-a)cos6+cos$„sin6D 

u  u  U  vjir^  0  0 


Cb12  b13] 


6B 

66Gr 
L  Grj 


+  C311  ai2] 


6A 


<5<i> 


(2-22) 

(2-23) 


(2-24) 


For  n  pairs  of  observations  one  obtains  the  nx2  matrix  A  and  the  nx2n  matrix  B. 


2.2.2  Case  b):  Horizontal  direction  and  time 


6B  must  be  eliminated  in  II.  For  that  purpose  equation  I  is  solved  for  5B: 


i3  ®ii  ^12 

6B  -  —  t- —  60  “  Tj —  ~  tt —  6$ 

°12  r  °12  12 


(2-25 


and  this  is  then  inserted  into  II: 


b216T+(b23_b22  ¥^)50Gr+(a2l'b22  F^)6A'b22  F^  ^+&236Z  =  ° 


(2-26 


Thereby  BQ  must  be  calculated  from  equation  I. 
Then  the  matrices  read: 


B  =  Cb  b  „  -  b 

j  21  23  u22  F^ 


(2-27] 


C-cosAQcosB0  -cosh0cos6-sinA0tanB0cos$0sinh0cos53 


a  a 

A  =  Ca  -b  T-ii  -b  a  1 

j  21  22  22  F^  23 


(2-28) 


C-cosh0cos6-sinA0tanB0cos<i>0sinh0cos6  -sinA^anB^sin^coshgCosa-cos^s^) 
-cosAQcosB^] 


h  h  13i 
b23  -b22 


+  Ca„  -b 


21  22  ^ 


"b22  F^-  a23]  6* 


(2-29 


For  n  pairs  of  observations  one  obtains  the  nx3  matrix  A  and  the  nx2n  matrix  B. 


2.2.3  Case  c):  Horizontal  direction,  vertical  direction  and  time 


Case  c)  is  represented  in  equations  (2-14)  to  (2-19). 


2.2.4  Case  d):  Horizontal  direction  and  vertical  direction 


S0Gr  must  be  eliminated  in  II.  For  that  purpose  equation  I  is  solved  for  60( 


60 


12 


Gr 


6B 


IP-  6A  -  ^  64> 
13  13 


and  then  this  is  inserted  into  II: 


b  a  a 

12  .  un./.  ll  u  12 


b216T+(b22"  B^  b23)6B+(a2r  B^  b23)6A"  B^  b236$+a236Z  =  0 


As  stated  above  a^  =  b13>  a  =  b23 


This  yields 


a  11  u 

21 '  23 


=  0! 


From  this  one  obtains  the  final  conditional  equation 


b2!5Tt<b22-  ^  b23M+a235!:  =  0 


Hence  it  follows  that  A  cannot  be  estimated.  The  reference  of  the  Greenwich 
meridian  for  0Qr  and  as  a  result  the  longitude  A  is  in  principle  chosen  arbi¬ 
trarily  so  that  A  is  not  definable  without  time  9  . 


The  matrices  read 


B:  ■  Cb21  b22  -  5^  V 


=  C-cosA.cosB.  sinA.sinB.  +  -  r  u  .  3 
oo  oo  cos$Qtanh0 


*J  *  C-  FT  b23  a233 


=  Ctan$0cosh0coth0cos6-sin6coth0  -cosA0cosBQ] 


For  n  pairs  of  observations  one  obtains  the  nx2  matrix  A  and  the  nx2n  matrix  B 


.VV.-.. V.VvV.V 


3. 


Numerical  studies 


3.1  The  adjustment  model 

As  shown  in  chapter  2  the  observations  and  unknowns  are  connected  with  each  other 
in  linearized  form  according  to  equation  (2-10)  or,  in  the  case  of  consistent  pa¬ 
rameters,  according  to  equation  (2-11).  This  leads  to  the  Gmb-HelmesU  model  of 
condition  equations  with  unknowns: 

Ax  +  Be  =  By  (2-11) 

In  the  realized  simulation  calculations  one  is  actually  not  interested  in  the  es¬ 
timation  of  the  unknowns  x,  but  only  in  the  accunncy  with  which  the  unknowns  can 
be  estimated  dependent  upon  the  actual  configuration  of  the  stars.  Thus  one  is 
only  interested  in  the  variance-covariance  matrix  of  the  unknowns.  To  determine 
this  matrix  the  coefficient  matrices  for  the  present  case  of  observations  are 
first  of  all  prepared:  the  parameter  matrix  A  and  the  condition  matrix  B.  Be¬ 
sides  this  the  variance-covariance  matrix  £,  i.e.  the  dispersion  matrix  D(y)=  I 
of  the  observation  vector  y,  has  to  be  stipulated  which  is  in  general  assumed  to 
be  a  diagonal  matrix. 

With  this  the  normal  equation  matrix  N  can  be  calculated: 

N  =  BZBt  (3-1) 

Finally  the  variance-covariance  matrix  Q  of  the  unknowns  follows  from  that  as 

XX 

Q.-DlV*-1  =  CAT(BlBTr1Ar 1  (3-2) 

The  square  roots  of  the  diagonal  elements  of  Q  are  now  the  accuracies  required 
with  which  the  unknowns  can  be  determined. 

3.2  Simulation  calculations 

Simulation  calculations  have  been  carried  out  for  the  four  different  possible  cases 
of  observations  of  stars.  Thereby  the  coordinates  right  ascension  a  and  declina¬ 
tion  6  of  fictitious  stars  were  determined  in  dependence  on  the  sidereal  time  in 
such  a  way  that  they  approximately  lay  in  defined  directions  as  seen  from  the  ob¬ 
servation  point.  In  accordance  with  this  the  simulation  calculations  were  carried 
out,  whereby  the  results  in  the  diagrams  are  linked  to  the  respective  star  configu- 
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rations  as  follows: 


the  stars  which  are  observed  here  are  positioned  approxi¬ 
mately  in  the  meridian  of  the  observation  point,  both  in 
the  South  as  well  as  in  the  North  (azimuth  s  0°  and  180° 
respectively) ; 

the  stars  are  positioned  approximately  in  the  first  verti¬ 
cal;  this  is  the  great  circle  through  the  directions  East 
and  West  (azimuth  s  90°  and  270°  respectively); 

u  t. 

the  stars  have  the  hour  angle  hs6  and  hs18  respectively, 
i.e.  considered  in  the  equatorial  frame  they  have  the  an¬ 
gular  distance  ±90°  from  the  meridian  of  the  observation 
point  measured  along  a  parallel  of  latitude; 

this  line  contains  stars  which  are  distributed  over  the 
whole  firmament:  meridian  and  first  vertical; 

the  stars  of  this  set  of  observations  are  located  partially 
in  the  meridian  and  partially  they  have  the  hour  angle  6h  and 
18h  respectively. 


In  the  illustrations  the  mean  errors  (or  standard  deviations)  of  the  unknowns  as¬ 
tronomical  longitude  A,  astronomical  latitude  $  and  orientation  unknown  of  the 
theodolite  I  (as  determinable  in  the  different  cases),  which  are  to  be  expected, 
are  drawn  in  dependence  upon  different  parameters.  For  the  four  different  cases 
of  observations  the  accuracies  in  the  determination  of  the  unknowns  are  given  first 
of  all  in  dependence  upon  the  number  of  the  observed  stars  and  subsequently  in  de¬ 
pendence  upon  the  accuracy  of  the  observations  (horizontal  and  vertical  direction 
or  time  measurement);  here  the  calculations  were  carried  out  with  ten  observed 
stars. 

For  the  calculation  of  the  accuracy  of  the  unknowns  with  the  free  parameter  "num¬ 
ber  of  stars"  the  accuracy  of  the  observations  has  been  supposed  as: 

horizontal  direction  :  a  =  1" 

T 

vertical  direction  :  a  ,=  1" 

time  measurement  :  oQ  =  0.1  sec 

The  accuracy  of  measuring  the  vertical  direction  a  ,  is  completed  by  the  accuracy 
of  the  determination  of  the  refraction,  which  can  be  calculated  in  dependence  upon 
the  vertical  direction  according  to  the  equations  in  Appendix  A.1.  These  results 
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differ,  however,  only  by  a  few  hundredths  of  a  second  of  arc  from  those,  which 
would  follow  when  one  neglects  the  inaccuracy  of  the  refraction.  In  some  illu¬ 
strations  the  ordinate  axis  is  drawn  as  a  broken  line.  This  means  that  the 
scale  in  the  upper  part  of  the  diagram  does  not  coincide  with  the  scale  in  the 
lower  part.  The  corresponding  values  are  explicitly  given. 

The  results  shown  in  the  different  illustrations  are  strictly  speaking  valid  only 
for  the  "observations"  which  have  been  supposed  here.  They  are  certainly  depen¬ 
dent  upon  the  constellation  of  the  stars  in  the  respective  group  of  the  observed 
stars.  The  tendency  of  the  results  will  be  correct  anyway. 

3.2.1  Case  a):  Vertical  direction  and  time 

Of  course  only  accuracies  for  longitude  and  latitude  could  be  calculated  here 
because  no  horizontal  directions  were  measured. 


matter  of  sear*'  n«Mr  of  stars 


Fig.  2  and  Fig.  3:  Accuracies  in  dependence  on  the  number  of  stars 
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F IguJiu  14,15  and  16:  Accuracies  in  dependence  upon  the  accuracy  of  the  time 

measurements  (10  stars,  oT  =  1") 


3.2.3  Case  c):  Horizontal  direction,  vertical  direction  and  time 


Figure  24  shows  that  in  this  case  an  inaccurate  time  measurement  has  hardly  any 
effect  on  the  determination  of  the  astronomical  latitude  (if  a  sufficient  number 
of  stars  have  been  observed). 
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Fig.  23,24  and  25:  Accuracies  in  dependence  upon  the  accuracy  of  the  time 

measurements  (10  stars,  cl,  =  o_,  =  1") 

3.2.4  Case  d):  Horizontal  direction  and  vertical  direction 

As  shown  in  2.2.4  the  astronomical  longitude  cannot  be  determined  without  measuring 
the  time;  thus  the  accuracy  of  the  longitude  cannot,  of  course,  be  specified  ei¬ 
ther. 
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Fig.  26  and  Fig.  27:  Accuracies  in  dependence  upon  the  number  of  stars 
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F-cg.  28  and  Fig.  29:  Accuracies  in  dependence  upon  the  accuracy  of  the  measured 

angles  (10  stars) 


3 . 3  Comi  ±r  programmes 


Two  computer  programmes  have  been  developed  in  FORTRAN  77,  which  are  implemented 
in  a  PDP11/23+  computer.  The  first  programme  APPROX  simulates  the  observation 
of  a  star.  After  the  input  of  approximate  values  for  the  astronomical  longitude 
and  latitude  of  the  observation  point,  the  orientation  unknown,  right  ascension 
and  declination  of  a  star  and  the  sidereal  time  the  programme  calculates  consis¬ 
tent  values  for  horizontal  and  vertical  direction  with  the  equations  (2-1)  to 
(2-3)  and  writes  the  complete  data  set  for  every  star  onto  a  data  file  APPVAL. 

By  choosing  the  respective  parameters  the  observation  directions  required  (see 
3.2)  can  be  computed. 


Ihe  second  programme  ASTRO  calculates  the  accuracies  with  which  the  unknowns  can 
be  estimated.  In  order  to  be  able  to  do  this  it  needs  two  data  files  as  the  in¬ 
put: 


f 


B 


I 

I 


s 

I 

* 


-  the  first  field  contains  the  chosen  data  sets  of  the  observed  stars  consisting 
of  the  observations  (according  to  case  a)  to  d))  and  the  coordinates  of  the 
star,  right  ascension  and  declination  (the  value  in  the  cases  a),  b)  and  d), 
which  is  not  an  "observation",  i.e.  horizontal  direction,  vertical  direction 
or  time,  completes  the  data  set  as  an  approximate  value); 

-  the  second  file  contains  three  approximate  values  for  the  unknowns  A,4>  and  Z. 

The  values  have  to  be  written  onto  the  file  in  the  format  deg(or  h)'(or  min)" 

(or  sec),  for  example  256  17  23.4  or  17  12  49.8.  In  Appendix  A. 2  the  first  pro¬ 
gramme  APPROX  and  the  second  programme  ASTROC  for  case  c)  are  listed. 


APPENDIX 


A.1  The  reduction  of  the  observed  vertical  direction  owing  to  refraction 
A. 1.1  The  reduction  formula 

The  vertical  direction  B  to  a  star  is  influenced  by  the  astronomical  refraction. 
Therefore  the  measured  vertical  direction  B1  must  be  corrected  by  the  value 
of  the  influence  of  refraction: 

B  =  B*-R  ( A— 1 ) 

R  is  the  refraction,  which  is  calculated  by  taking  as  a  basis  an  atmospheric 
model  and  using  temperature  and  atmospheric  pressure  measurements.  The  first  step 
for  calculation  is  the  determination  of  the  normal  refraction  RQ  which  is  valid 
for  the  following  conditions  at  the  observation  point: 

atmospheric  pressure  pQ  =  760Torr  =  1013,25  hpa 
temperature  tQ  =  0°C 

vapour  pressure  eQ  =  6Torr 

According  to  K.  Ramsayer  CHandbuch  der  Vermessungskunde,  Band  I la :  Geodatische 
Astronomie,  1968,  p.  115ffl  the  normal  refraction  can  be  calculated  with  high  pre¬ 
cision  by  the  approximation 

RQ  =  60", 1012  •  cotB  -  0", 06483  •  cot3B  (B>20°)  (A-2) 

This  equation  is  based  on  refraction  indices  of  the  atmosphere  which  are  calculated 
for  different  altitudes  with  a  mean  distribution  of  atmospheric  pressure  and  tem¬ 
perature,  which  by  the  way  is  in  good  agreement  with  the  U.S.  Standard  Atmosphere 
1962. 


In  the  second  step  there  follows  the  reduction  to  another  atmospheric  pressure  p 
and  to  another  temperature  t  at  the  observation  point.  The  equation  for  this 
reads  as  follows: 


and  with  the  abbreviations 


pChpa] 

l0l3,25hpa 


R  =  RQ  •  G  •  K 


„  273,15 

iP'.i 


The  reduction  to  another  vapour  pressure  e  is  calculated  as  follows 


dR  =  -0" ,011 (e-6)  •  K  •  cotB 

or  in  consideration  of  equation  (A-2) 

dR  =  -0,00018  •  R  •  K  •  (e-6) 

This  reduction  can  generally  be  neglected. 


A. 1.2  The  accuracy  of  the  reduced  vertical  direction 


The  accuracy  of  the  vertical  direction  depends  firstly  on  the  measuring  accuracy 
and  secondly  on  the  accurcy  of  the  calculation  of  the  refraction.  This  problem 
will  now  be  studied. 

A. 1.2.1  Deviations- from  the  normal  state  of  the  atmosphere 


According  to  K.  Ramsayer  the  equation  of  refraction  is  valid: 


l  n 


.  1 


R0  =  p(n  -1 ) •cotB-cotB*j*  /  h«d(-^-)  +  ^*p*(n0-1  )2*cot3B 


**  l  n 

-cot3B*~-*  f  h -d(— P-) 
a  n 

no 


o  =  (180*  3600  >/tt  ["] 


refraction  index  at  the  observation  point 

radius  of  the  osculating  sphere  at  the  observation  point  ("earth  radius") 
refraction  index  of  a  stratum  of  the  atmosphere  in  the  altitude  h. 


n 


o 


a 


In  the  case  of  the  mean  state  of  the  atmosphere  equation  (A-7)  yields  with 
pQ  =  760,3Torr,  tQ  =  9,4°C,  eQ  =  4,8Torr  at  the  observation  point  the  refrac¬ 
tional  value  as  follows: 

RQ  =  58",282*cotB-0",076*cotB+0",0082*cot3B-0",0762‘COt3B  .  (A-8) 

The  comparison  of  equation  (A-7)  and  (A-8)  proves  that  the  refraction  primarily 
depends  on  the  refraction  index  at  the  observation  point  and  that  it  is  to  a 
large  extent  independent  of  the  change  of  the  refraction  index  which  occurs  with 
the  change  of  the  altitude. 

Supposing  a  horizontal  and  plane  stratification  of  the  atmosphere  (a-»°)  equa¬ 
tion  (A-7)  yields 

Rq  *  p(nQ-1 ) *cotB+  ^(n0-1)2«cot3B  (A-9) 

Equation  (A-9)  shows  that  in  the  case  of  this  assumption  (a-*=°)  the  refraction  is 
only  dependent  on  the  vertical  direction  and  the  refraction  index  at  the  obser¬ 
vation  point  and  is  completely  independent  of  the  state  of  the  atmosphere  above 
it.  The  fact  that  the  state  of  the  atmosphere  does  have  an  effect  on  the  refrac¬ 
tion  all  the  same  (in  accordance  with  the  integral  terms  in  equation  (A-7))  is 
the  consequence  of  the  curvature  of  the  optical  strata  of  the  atmosphere.  The 
integral  terms  have  the  following  value  for  a  mean  state  of  the  atmosphere 

0" ,076*cotB  and  0" ,0762-cot3B  , 

and  for  B>30°  they  reach  an  order  of  magnitude  of  only  a  few  tenths  of  a  second 
of  arc.  Hence  it  follows  with  the  assumption  that  the  optical  strata  are  concen¬ 
tric  spheres  that  even  considerable  deviations  from  the  normal  state  of  the  at¬ 
mosphere  (for  example  caused  by  a  temperature  inversion  in  the  lower  strata) 
can  change  the  refraction  at  the  most  by  0",1  to  G",2. 


So  an  estimation  for  the  loss  of  accuracy  of  the  refraction  because  of  the  de¬ 
viation  from  the  normal  state  of  the  atmosphere  might  be 

dRt  *  0",05  •  cotB  (B>20°)  (A-10) 

or  for  the  expected  mean  error  of  the  determination  of  refraction  because  of 
the  deviation  from  the  normal  state  of  the  atmosphere 

m  =  ±0",05  •  cotB  (B>20°)  (A-lla) 

This  is  equal  to  the  standard  deviation: 

ot  =  0",05  •  cotB  (B>20°)  (A-llb) 

A. 1.2. 2  Inclination  of  the  strata  and  deviation  from  the  spherical  form 

Most  of  the  refraction  theories  postulate  that  the  optical  strata,  i.e.  the 
strata  with  the  same  refraction  index  of  the  atmosphere,  are  concentric  spheres 
which  are  perpendicular  to  the  plumb-line  at  the  observation  point.  In  reality 
these  assumptions  are  not  exactly  complied  with.  One  has  to  rather  reckon  with 
an  inclination  of  the  strata  and  a  deviation  from  their  spherical  form. 

An  estimation  of  this  influence  yields 

AR  *  0"  ,084  •  hCkm]  •  ijl  C^]  •  cosec2B  (A-12) 

The  altitude  of  adjustment  h  specifies  up  to  which  altitude  one  has  to  reckon 
with  an  inclination  of  the  strata,  is  the  horizontal  temperature  gradient. 

From  series  of  observations  which  Harzer  C'Berechnung  der  Ablenkungen  der  Licht- 
strahlen  in  der  Atmosphere  der  Erde  auf  rein  meteorologisch-physikal ischer  Grund- 
lage",  Publ .  der  Sternwarte  in  Kiel  XITI,  Germany,  1922-24]  carried  out, the  influ¬ 
ence  of  the  deviation  of  the  mean  real  optical  surfaces  from  a  concentric  spheri¬ 
cal  form  results  in  0",01  for  B  =  70°  and  0",Q6  for  B  =  30°  respectively. 

A  valuation  for  this  error  might  therefore  be 
o2  =  0",015*cosec2B 


(A- 13 ) 
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A. 1.2. 3  The  influence  of  the  real  water  vapour  pressure 

The  formula  for  the  normal  refraction  uses  a  water  vapour  pressure  of  eQ  =  6Torr. 
According  to  K.  Ramsayer  the  influence  of  the  real  vapour  pressure  e  is 

AR  =  0“011(eo-e)*cotB  .  (A- 1 4) 

e  oscillates  seasonably  between  4Torr  in  January  and  12Torr  in  August  (in  Germany!). 

The  influence  thereby  becomes 

I AR I  *  0" ,07  •  cotB  (A- 15) 

Generally  it  can  be  neglected.  The  expected  mean  error  or  the  standard  devi¬ 
ation  could  be  stipulated  to 

o3  *  0" ,02  •  cotB  (A- 16) 

Hereby  one  supposed  that  in  general  e  does  not  vary  much  from  eQ. 

A. 1.2. 4  The  influence  of  the  wave-length  of  the  starlight 

The  refraction  is  dependent  on  the  wave-length  of  the  light  emitted  by  the  stars. 

The  change  of  the  refraction  varies  between  +6" ,05*cotB  (white  fixed  stars)  and 
-0",23«cotB  (red  fixed  stars).  Thus  the  colour  correction  amounts  to  some  tenths 
of  a  second  of  arc  for  B<45°  in  the  case  of  reddish  yellow  to  red  stars,  and 
should,  therefore,  be  allowed  for  when  calculating  the  refraction.  An  estimation 
for  the  mean  error  is,  therefore,  not  made. 

A. 1.2. 5  The  influence  of  errors  in  the  measurements 


Finally,  the  influence  of  an  error  in  measuring  atmospheric 
ture  must  also  be  considered.  Starting  from  equation  (A-3) 
lowing  when  assuming  the  values  R  *  60"*cotB,  p  =  760Torr,T 


pressure  and  tempera- 
one  obtains  the  fol- 
=  283°K: 


dR  *  0" ,079*cotB*  dp  -  0",212*cotB*dT 


(A— 1 7) 


The  change  to  the  standard  deviation  yields  with  op  =  0,1Torr  and  aT  =  0,1°C 
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So  the  influence  of  an  error  in  pressure  and  temperature  can  in  general  be  ne¬ 
glected. 

A. 1.2. 6  Summary  of  the  different  errors 

Summing  up  the  proceeding  influences  on  refraction  this  yields 


/  2  2  2  2 
°R  =  /ai+a2+°3+a5 


=  /  C(0" ,05)2+(0" ,02)2+(0" ,023)2]cot2B+(0"  ,015)2cosec4B 
and  finally 

a  =  /  (0\06-cotB)2+(0",015-cosec2B)2  . 


( A— 1 9 ) 
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This  error  has  to  be  added  to  the  observational  error  oQ  from  which  the  standard 
deviation  of  the  corrected  vertical  direction  B  results: 

°B  *  Jal  (fl-20> 

A. 1.2. 7  Comparison  with  series  of  .observations 

Im  comparison  with  equation  (A-19)  the  results  of  series  of  observations  of  meri¬ 
dian  vertical  directions  lasting  one  and  two  years  carried  out  by  J.  Bauschinger 
(1898)  and  L.  Courvoisier  (1904)  are  presented.  Averaging  the  differences  between 
calculated  and  observed  vertical  directions  one  obtains 

Bauschinger  :  a  =  /  (0"  ,31 )2+(0" ,0035*cosec2B)2  (A-21) 

Courvoisier  :  a  *  / (0" ,23)2+(0",032*cosec2B)2  (A-22) 


The  first  terms  in  the  square  root  denote  the  influence  of  the  observation  errors. 
The  second  terms  can  be  interpreted  as  the  mean  influence  of  the  inclination  of 
the  optical  strata  and  represent  here  the  real  error  of  the  determination  of  the 
refraction.  So  the  value  in  equation  (A-13)  assumed  for  the  influence  of  the  in- 


clination  of  the  optical  strata  lies  between  those  values  which  were  deducted 
from  observations  in  equations  (A-21)  and  (A-22).  The  parts  of  the  errors  as 
in  equations  (A-11)  and  (A~1 6)  do  not  appear  here  since  convenient  conditions 
were  probably  present  for  the  observations.  The  part  of  the  error  from  equa¬ 
tion  (A-18)  could  be  eliminated  by  observing  diametric  stars  lying  more  or  less 
in  the  same  horizontal  distance  and  by  averaging  these  observations. 
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A. 2  Programme  listings 


A. 2.1  Programme  APPROX 


PROGRAM  APPROX 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

REAL  MM 

CHARACTERS  DUMMY 
PI=4.*ATAN(  1 . ) 

OPEN(  1  ,FILE='APPVAL',ST  ATUS='NEW' ) 

1  WRITE(5,100) 

100  FORMAT  (5X, 'astronomical  longitude  (deg  "  ")') 

READ(5,*)G,BM.  BS 

CALL  OUT(G,BM,BS,AL,1 ) 

WRITE(5,101) 

101  FORMAT!  5X, 'astronomical  latitude  (deg  "  ")') 

READ(5,*)G,BM,  BS 

CALL  OUT(G,BM.BS,AB,2) 

WRITE(5, 102 ) 

102  FORMAT (5X, 'orientation  unknown  (deg  "  " )') 

READ(5,*)G,BM,BS 

CALL  OUT(G.BM,BS,OU,3) 

WRITE(5,103) 

103  FORMAT(5X, 'right  ascension  (h  min  sec)') 

READ(5,*)H,MM,SS 

CALL  OUT(  H.MM.SS,  RA,4) 

WRITE(5,104) 

104  FORMAT ( 5X,  'declination  (deg  "  ")') 

READ(5,*)G,BM,BS 

CALL  OUT(G,BM.BS,D,5) 

WRITE(5, 105) 

105  FORMAT!  5X, 'Greenwich  sidereal  time  (CAST)  (h  min  sec)') 
READ(5,*)H,MM,SS 

CALL  OUT(  H,MM,SS,TH,6) 

SINB=DCOS(  AB  )*DCOS(  TH+AL-RA  )*DCOS(  D)+DSIN(  AB)*DSIN(  D) 
B=DASIN(SINB) 

SI NAZ=-DSIN( TH+AL-RA  )*DCOS(  D )/DCOS(  B ) 

AZ1=DASIN(SINAZ) 

COSAZ=( DSIN(  AB  )*DCOS(  TH+AL-RA  )*DCOS(  D  )-DCOS(  AB  )*DSIN(  D ) ) 
+  /DCOS(B) 

AZ2=DACOS( COSAZ ) 

IF(SINAZ.GE.O. .  AND.  COSAZ .  CE.0. )  AZ=AZ2 
IF(SINAZ.CE.O.. AND. COSAZ. LT.0. )  AZ=AZ2 
IF(SINAZ.LT.0.  .AND. COSAZ. LE.0. )  AZ=2.*PI-AZ2 
IF(SINAZ.LT.O.  .AND. COSAZ. CT.0. )  AZ=2.*PI-AZ2 
T=AZ-OU 

I F (  T .  LT.0. )  T=T+2.*PI 
Output  to  the  terminal 

WRITE(5, '(///, 5X, "longitude:", F12. 6,"  deg")')  AL*180./PI 
WRITE(5,'(5X. "latitude:", F12. 6,"  deg")')  AB*180./PI 


uuuuu 
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WRITE(5,'(5X, "orientation  unknown:", F12. 6,"  deg")') 

+  OU*180./PI 

WRITE(5,'{5X, "right  ascension:"F12.6,"  h" )' )RA*180. /15. /PI 
WRITE(5,'(5Xy "declination :",F12 .6,"  deg")')  D*180./PI 
WRITE(5,'(5X, "sidereal  time:",F12.6,"  h")')  TH*180. /15. /PI 
WRITE(5,'(5X, "vertical  direction:", F12. 6,"  deg" )' )B*180. /PI 
WRITE(5,'(5X, "horizontal  direction:", F12. 6,"  deg")') 

+  T*180./PI 

WRITE(5,'(5X, "azimuth:". F12. 6."  deg")')  AZ*180./PI 

WRITE(5,'(5X, "sin (azimuth): ",F12.6)')  SINAZ 

WRITE(5,'(5X,"cos(azimuth)  :",F12.6)')  COSAZ 

B=B*180./PI 

T=T*180./PI 

AZ=AZ*180./PI 

IC=INT(  B ) 

I F ( B . LT . 0 . )  IG=-IG 
BM=ABS(B-IG)*60. 

IBM=INT(BM) 

BS=(BM-IBM)*60. 

WRITE!  1,1 06 )  I G ,  IBM. BS 

106  FORMAT! 5X, 'vertical  direction: ',14, 'deg',13, "", F5 . 1 ,"" ) 
IG=INT(T) 

BM=(T-IG)*60. 

IBM=INT(BM) 

BS=  ( BM- 1 BM )  *60 . 

WRITE(1,107)IG,IBM,BS 

107  FORMAT(5X, 'horizontal  direction: ',14, 'deg',  13, "",F5. 1 ) 
IG=INT(  AZ) 

BM=(AZ-IG)*60. 

IBM=INT(  BM) 

BS=  ( BM-  IBM)  *60 . 

WRITE!  1,1  OS  JIG,  IBM,  BS 

108  FORMAT(5X, 'azimuth: '.14,'deg',  I3,"",F5.1,"") 

WR I TE(1,'(5X, "sin! azimuth)  :",F12.6,5X, "cos! azimuth) :", 

+  FI 0 . 6 ) ' )  SINAZ, COSAZ 

WRITE(5,'(5X, "Should  a  new  data  set  be  calculated  (Y/NU  ? 

+  ")') 

READ(5,'(  A1 )' )  DUMMY 

IF(DUMMY.EQ.'Y'.OR.DUMMY.EQ.'y'.OR  DUMMY. EQ. '1 ' )  GOTO 

STOP 

END 
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200 

10 

201 

20 
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40 
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60 


SUBROUTINE  OUT(G,BM,  BS, RET,N ) 

DOUBLE  PRECISION  G,  BM,  BS,  RET, PI 
Output  to  the  data  file  APPV/AL 
PI=4.0*ATAN(  1.0) 

IG=INT(  G) 

IBM=INT(BM) 

IF(N.GT.I)  GOTO  10 
WRITE!  1 ,200)  IG,  IBM,  BS 

FORMAT! / / / , 5X , 'astronomical  longitude: ',  14, 'deg1, 13,"", 

+  F5  1  1111 ) 

RET=(G+BM/60.+BS/3600.  )*PI/180. 

GOTO  60 

IF(N.GT.2)  GOTO  20 
WRITE!  1 ,201 ) IG,  IBM, BS 

FORMAT!  5X, 'astronomical  latitude:',  14, 'deg',  13, "",F5. 1 ) 
RET=l G+BM/60.+BS/3600.  )*PI/180. 

GOTO  60 

IF(N.GT.3)  GOTO  30 
WRITE!  1,202 )  IG,  IBM,  BS 

FORMAT!5X, 'orientation  unknown: ',  14, 'deg*,  13,"", F5. 1 ,"" ) 
RET=(G+BM/60.+BS/3600.  )*PI/180. 

GOTO  60 

IF(N.GT.4)  GOTO  40 
WRITE!  1 ,203)  I G,  IBM,  BS 

FORMATI5X,  'right  ascension:',  14,'h',  l3,'min',F5. 1,'sec') 
RET=(G+BM/60.+BS/3600. ) *PI / 180. *15 . 

GOTO  60 

IF! N . GT.5 )  GOTO  50 
WRITE!  1 ,204)  I G,  IBM, BS 

FORMATI5X,  'declination:  ',14,  'deg',  13,  "",F5.1,"") 

RET=( G+BM/60.+BS/3600.  )*PI/180. 

GOTO  60 
CONTINUE 

WRITE!  1,205  JIG,  IBM, BS 

FORMAT ( 5X, 'Greenwich  sidereal  time:',  14,'h', 1 3, 'min', F5. 1 , 'sec' ) 
RET=(C+BM/60.+BS/3600.  )*15.*PI/180. 

RETURN 


PROGRAM  ASTROC 
REAL  MM 

PARAMETER  (KK=10) 

DIMENSION  SIGMA(3*KK),U(3),DISP(3,3) 

DIMENSION  Y(KK,3),SIGY(3*KK,3*KK),A(2*KK,3), 
B(2*KK,3*KK),BT(3*KK,2*KK),AT(3,2*KK),RN(2*KK,2*KK), 
ATN(3,2*KK),HILF(2*KK,3*KK),NUM(2*KK),D(2*KK).DELTA(  KK) 
ALPHA! KK) 

CHARACTER  FRACE*1 
DATA  RN/400*0. DO/ 


CHARACTER  *25  FNAME1 , FNAME2 
LOGICAL  LOG 


H=0. 

G=0. 

MM=0. 

BM=0. 

SS=0. 

BS=0. 

PI=4.*AT AN!  1 .  ) 

WRITE(5.'(//,5X,"name  of  the  data  file?")') 

READ(5,'(  A25)')FNAME1 

WR  ITE(5.'( //,5X,"name  of  the  file  with  the  approximation  values" 
"  for  the  unknowns  ?")') 

READ! 5, '{ A25)')FNAME2 

WRITE(5,'(//,5X, "accuracy  of  horizontal  directions  (")  ?")') 

READ(5,*)BS 

GHR=BS 

CALL  DEZIG!  G,BM,BS,SIGMA(  1 ) ) 

CALL  RADI  (SIGMA!  1 ) , P I ) 

WRITE(5,'(  /  /,5X,  "accuracy  of  vertical  directions  (")  ?")') 

READ(5,*)BS 

GHD=BS 

WR I TE  ( 5, '(//,  5X,  "consideration  of  the  accuracy  of  refraction", 

"  (Y/N)  ?")') 

READ(5,'(  A1 )'  )FRAGE 

WRITE(5,'( //,5X, "accuracy  of  time  observations  (sec)  ?")') 

READ(5,*)SS 

GZ=SS 

CALL  DEZIH(H,MM.SS,SIGMA(3)) 

CALL  RADI  (SIGMA!  3 ) ,  PI ) 
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OPEN(  1  ,FILE=FNAME1  ,  ACCESS='SEQUENTI  AL',STATUS='OLD' ) 
OPEN(2,FILE=FNAME2,  ACCESS='SEQUENTI  AL',ST  ATUS-'OLD' ) 
0PEN(3,FILE='0UTC.  DAT',ST  ATUS=,NEWI ) 


READ(2,*)G,BM,BS 

IG=INT(G) 

IBM=INT(BM) 

CALL  DEZIG(G,BM,BS,U(1 ) ) 
CALL  RADI(U(1),PI) 
READ(2,*)G,BM,BS 
IC=INT(G) 

IBM=INT(BM) 

CALL  DEZIG(G,  BM,  BS. U(2 ) ) 
CALL  RADI  ( U(2 )  ,PI ) 
READ(2.*)G.BM.BS 
IC=INT(C) 

IBM=INT(  BM) 

CALL  DEZIG(G,BM,BS,U(3) ) 
CALL  RADI  ( U(3),PI ) 


WRITE(3,700) 

700  FORMAT!///, 5X, 'observations', /,5x, 'turn:  horizontal  ', 
+  'direction  ', '(deg-"-")', /.11x, 'vertical  direction  ', 

+  '(deg-"-")', / ,  1 1 X , ‘sidereal  time  (GAST)  (h-min-sec)', 

+  /.  11X, 'and', / ,  1 1 X , 'right  ascension  (h-min-sec)', /,11X, 

+  'declination  (deg-"-")') 

K=1 

25  READ(1  ,  *,END=20) 

K=K+1 
GOTO  25 
20  K=K/5 

REWIND  1 


DO  10  J=1 ,  K 
WRITE(3. '(/)') 

DO  30  1=1,2 
READ(  1,*)G,BM,BS 
1  G=  I  NT  ( G ) 

IBM=INT(  BM) 
WRITE(3,800)  IG,  IBM,  BS 


IF!  I.EQ.2)  THEN 

BS=SQRT(GHD*GHD+(0.06/TAN(  Y(  J,  I ) )  )**2.+ 

+  (0.015/(S1N(Y(J, I ) ) )**2. )**2. ) 

IF(FRAGE.  EQ.'N'.OR  .FRAGE.  EQ.'n'  .OR  .FRAGE.EQ.  'O'  )THEN 

BS=GHD 

ENDIF 

CALL  DEZIG(0.  ,0. ,  BS,  SIGMA!  3*!  J-1  )+2 ) ) 

CALL  RADI(SIGMA(3*(J-1)+2),PI) 

ELSE 

SIGMA(3*(  J-1  )+1  )=SIGMA(  1 ) 

ENDIF 

3  CONTINUE 

READ!  1  ,*)H,MM,SS 
IH=INT(H) 

IMM=INT(MM) 

WRITE(3,'(5X,  14,  "h“,  13, "min", F5 . 1 ,  "sec" )' )  IH,  IMM.SS 
CALL  DEZIH(H,MM,SS,  Y(  J,3) ) 

CALL  RADI  ( Y(  J,3)  ,PI ) 

SICMA(3*(  J-1  )+3)=SIGMA(3) 

READ(  1  ,*)H.MM,SS 
IH=INT(H) 

IMM=INT(MM) 

WRITE(3,  '(5X,  14, “h1 1 , 13,  "min", F5. 1 ,  "sec" )' )  IH,  IMM.SS 
CALL  DEZ!H(  H.MM,SS,  ALPHA!  J  ) ) 

CALL  RADKALPHA(J).PI) 

READ!  1  ,*)G,BM,BS 
IG=INT(G) 

IBM=INT(  BM) 

WRITE(3,800)IG,IBM,BS 

CALL  DEZIG(G.BM,BS, DELTA!  J ) ) 

CALL  RADI  ( DELTA!  J  )  ,PI ) 

3  CONTINUE 

WRITE!  3, 801  )GHR 

I  FORMAT! /,5X, 'accuracy  of  horizontal  directions:', F6. 2, 

+ 

WRITE(3,802)GHD 

?  FORMAT(5X,'accuracy  of  vertical  directions  (observations):1 
+  ,F6.2,"" ) 

WRITE(3,803)GZ 

}  FORMAT(5X,'accuracy  of  time  observations: F6. 2, 'sec' ) 


DO  11  1=1 ,3*KK 
DO  11  J  =  1 ,3*KK 
SIGY(  I,  J  )=0. 

DO  12  1=1 ,2*KK 
DO  12  J=1 ,3 
A(  I,  J  )=0. 

DO  13  l  =  1,2*KK 
DO  13  J=1 ,3*KK 
B(  I,  J  )=0. 
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DO  14  1=1. 3*K 

14  WRITE(3.*)SIGMA(I)*180./PI*3600. 

CALL  SIGMAS(SIGY,KK,K, SIGMA) 

SIGMA  =  vector  of  the  standard  deviations  of  the  observations  (in  (rad)) 
SI  GY  =  matrix  of  the  variances  of  the  observations 
FAKTOR=SIGY(1.1) 

DO  32  1=1. 3*K 
DO  32  J=1.3*K 

SIGY(  I ,  J  )=S  I  GY  ( I ,  J  l/FAKTOR 
32  CONTINUE 

CALL  AMAT( A, Y.U.KK.K. ALPHA. DELTA) 

DO  40  1=1. 2*K 

WRITE(3,300)(A(I.M).M=1.3) 

100  FORMAT!  5X.3F  18. 12 ) 

40  CONTINUE 

CALL  BMAT(B. Y.U.KK.K, ALPHA, DELTA) 

1=0 

DO  41  J=1 ,2*K-1 ,2 
WRITE(3,300) ( B(  J,  l*3+M),M=1 ,3) 

WR  ITE(  3,300 )( B(J+1 .  I*3+M)  .M=1 ,3) 

WRITE(3, '(/)') 

41  1=1+1 


CALL  M ATMU L  ( B ,  2*  K  K ,  2*  K ,  3*  K  K ,  3*  K ,  S I G  Y ,  3*  K  K ,  3*  K ,  3*  K  K ,  3*  K , 

+  HILF,2*KK.2*K,3*KK,3*K ) 

DO  50  1=1. 2*K 

DO  50  J=1 ,3*K 

50  BT(J,  I  )=B(  I.J  ) 

CALL  MATMUL(  HILF.  2*KK,2*K,  3*KK,3*K,  BT,3*KK,3*K,2*KK,2*K,  RN, 
+  2*KK,2*K,2*KK,2*K ) 
normal  equations  matrix 

WRITE(3,'(5X. "matrix  B*SIGY*BT :")') 

DO  51  1=1. 2*K 

51  WRITE(  3,401  )(RN(I.J),J=1,2*K) 

CALL  INVER2(  RN,2*K,  NUM,  LOG,  D,  2*KK ) 

WRITE(3,'(5x. "determinant  of  (  B*SICY*BT) 

+  3X.E18. 10)'  )D(  1 ) 

WRITE(3,'(5X. "matrix  (B*SICY*BT)  -1:")') 

DO  52  1=1, 2*K 

52  WRITE(  3,401  )(RN(I,J),J=1,2*K) 

01  FORMAT(8E15.6) 

DO  60  1=1. 2*K 

DO  60  J=1 ,3 

60  AT(  J,  I  )=A(  I ,  J ) 

CALL  MATMUL( AT ,3,3,2*KK,2*K,RN,2*KK,2*K,2*KK,2*K,ATN,3,3, 

+  2*KK,2*K) 

CALL  MATMUL( ATN,3,3,2*KK,2*K, A,2*KK,2*K,3/3, DISP.3,3,3,3) 

CALL  INVER2(  DISP,  3.  NUM,  LOG,  D,3) 

WRITE(3.'(5x, "determinant  of  ( AT*(  ( B*SIGY*BT )  '-1  )*A ) : " , 

+  3X.E18.10)’)D(1) 
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variance-covariance-matrix  of  the  unknowns 
DO  70  1=1.3 

70  WRITE(3,500)(DISP(I.J)*FAKTOR,J=1,3) 

500  FORMAT!  3F22. 16) 

DO  90  1=1.3 

DISP(  1.1  )=SQRT(  DISP!  1,1 )  )*180. /PI*SQRT(FAKTOR) 
DISP(  I,  I  )=DISP(  I ,  I  )*3600. 

standard  deviations  of  the  unknowns  (in  seconds  of  arc) 
90  WRITE(3,600)  DISP(IJ) 

600  FORMAT  (5X,  FI  2. 4,'  "') 

STOP 

END 


SUBROUTINE  DEZIG! G. BM, BS, RET )  . 

DEZIG  transforms  the  input  ( degree, minute  of  arc. second  of  arc)  j) 

into  a  decimal  value  (degree)  j  » 

RET=G+BM/60.+BS/3600.  > 

RETURN 

END  u 


SUBROUTINE  DEZIH(  H,MM,SS,  RET ) 

DEZIH  transforms  the  input  (hour, minute  of  arc, second  of  arc) 
into  a  decimal  value  (degree) 

REAL  MM 

RET=( H+MM/60.+SS/3600.  )*15. 

RETURN 

END 


SUBROUTINE  RADI(RET.PI) 

RADI  transforms  from  (degree)  to  (rad) 
RET=RET*PI/180. 

RETURN 

END 
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SUBROUTINE  SICMASISIGY.  KK.  K, SIGMA ) 

SIGMAS  constructs  the  variance-covariance-matrix  of  the  observations 
DIMENSION  SI  GY(  3*KK,  3*KK ),  SIGMA  (3*KK ) 

DO  10  1=1, 3*  K 

SICY(  I ,  I  )=SIGMA(  I  )*SIGMA{  I ) 
l  CONTINUE 
RETURN 
END 


SUBROUTINE  AMAT!  A,  Y,U,  KK,  K,  ALPHA, DELTA) 

AMAT  constructs  the  matrix  A 

DIMENSION  A{2*KK,3),  Y(  KK,3 ), U( 3)  .ALPHA!  KK ), DELTA!  KK ) 

J=1 

DO  10  1=1 ,2*K— 1 ,2 

A  ( 1 , 1  )=-COS(  U(2)  )*SIN(Y!(l+1)/2,3)+U(1 ) -ALPHA  ( J )  )*COS(  DELTA!  J ) ) 
A(  1 .2  )=-SIN(  U(2 )  )*COS(  Y(  ( 1+1 ) /2,3 )+U(  1 ) -ALPHA!  J  )  )*COS(  DELTA!  J  ) ) 
+  +COS!  U|2)  )*SIN(  DELTA!  J ) ) 

J=J+1 

CONTINUE 

J=1 

DO  20  1  =2 , 2*  K ,  2 

A(  L  1  )=-COS(  Y(  1/2,3 )+U  ( 1  )-ALPHA(  J )  )*COS(  DELTA!  J  ) ) 

A(  l.3)=-COS(U(3)+Y (1/2,1 )  )*COS(  Y(  1/2,2) ) 

J=J+1 

CONTINUE 

RETURN 

END 


SUBROUTINE  BMAT(  B,Y,U,  KK,  K,  ALPHA, DELTA ) 

BMAT  constructs  the  matrix  B 

DIMENSION  B(2*KK,  3*KK).Y(KK.3),U(  3),  ALPHA!  KK),  DELTA!  KK) 
1=0 

DO  10  J=1 ,2*K-1 ,2 
B ( J , 1*3+1 )=0. 

BlJ,  1*3+2 )=-COS(Y(  1+1,2) ) 

B(  J.  1*3+3 )=-COS( U (2 )  )*SIN!  Y(  1  +  1 ,3)+U(  1 ) -ALPHA!  1+1 ) ) 

•  *COS(  DELTA!  1+1 )) 

1  =  1+1 

CONTINUE 

1=0 

DO  20  J=2,2*K,2 

B(  J.  1*3+1  )=-COS(U(3)+Y(  1+1, 1))*COS(Y(  1+1,2)) 

B(  J.  1*3+2  )=SIN!U(3)+Y(  1+1. 1))*SIN(Y(  1+1,2)) 

B ( J ,  1*3+3 )=-COS(Y(  1+1.3 )+U ( 1 J-ALPHA!  1+1 )  )*COS(  DELTA!  1  +  1 ) ) 

1  =  1+1 

CONTINUE 

RETURN 

END 
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SUBROUTINE  MATMULl  RM1 ,  NDZ1 .  NZ1 ,  NDS1 ,  NS1 ,  RM2,  NDZ2,  NZ2, 
+  NDS2,NS2,RM3,NDZ3,NZ3.NDS3,NS3) 

MATMUL  executes  the  product  of  two  matrices 

DIMENSION  RM1  ( NDZ1 ,  NDS1) , RM2( NDZ2, NDS2 ) , RM3I NDZ3,  NDS3) 
DO  10  1=1, NZ3 
DO  10  J=1.NS3 
RM3(  I , J )=0. 

DO  10  11=1, NS1 

RM3(  I.  J  )=RM3(  I,  J  )+RM1  (1,11  )*RM2(  1 1,  J ) 

3  CONTINUE 
RETURN 
END 


SUBROUTINE  INVER2  ( B,N,NUM,LOC,D,NA) 
INVER2  executes  the  matrix  inversion 
DIMENSION  B(NA.NA),NUM(NA),D(NA) 


CALL  NORMEN  (B,N,N,ZB.SB,QB, AB,NA) 


DO  1  1=1,  N 
1  NUM(  I  )=0 
DET=1 . 

DO  2  L=1.N 
DO  3  M=1,N 

IF  (NUM(M).NE.O)  GOTO  5 
DIM)  =  ABS  (B(M,M) ) 

GOTO  3 

5  DIM)  =  0. 

3  CONTINUE 

II  =.1 

DO  6  1=2, N 

IF  (DUD. LT. Dll))  11  =  1 

6  CONTINUE 

IF  (DIID.GE.1.E-30)  GOTO  31 
NUM  (ID  =  L 
GOTO  2 

31  NUM  (ID  =  L 
DO  10  J  =  1 ,  N 
IF  (J.EQ.I1)  GOTO  10 
DO  13  1=1. N 
IF  (I.NE.I1)  GOTO  20 
DUD  =  BU1.J)/BU1.ID 
GOTO  13 

20  Dll)  =  B(  I,  J)-Bll.M  )*B(I1,J)/B(I1 . 11) 

13  CONTINUE 
DO  14  1=1, N 

14  BII.J)  =  DU) 


m 


CONTINUE 
H  =  B  ( 1 1 , 1 1 ) 

DET  =  DET*H 
00  15  1=1, N 
IF  (I.NE.I1)  GOTO  17 
B(  I,  II )  =  1/H 
GOTO  15 

B(  I,  II )  =  -B(  I ,  II  )/H 

CONTINUE 

CONTINUE 

K=0 

DO  18  1=1, N 
K=K+NUM(  I ) 

KV  =  CN*(N+1 )  )/2 
LOG  =  .FALSE. 

IF  (K.EQ.KV)  LOG=. TRUE. 


CALL  NORMEN  (B,N,N,ZB1,SB1,QB1,AB1,NA) 

IF  (N.GE.6)  D(6)=AC*AB1 

IF  (N.GE.5)  D(5)=QB*QB1 

IF  (N.GE.4)  D(4)=SB*SB1 

IF  (N.CE.3)  0{3)=ZB*ZB1 

IF  (N.GE.2)  D(2 )=1  /DET 

D  ( 1  )=DET 

RETURN 

END 


SUBROUTINE  NORMEN  ( A,N,M,Z,S,Q.AC,NA) 
DIMENSION  A(NA,NA) 


Z=0 

DO  1  J=1,N 
ZS=0 

DO  2  K=1,M 
ZS=ZS+ABS  ( A(  J,  K ) ) 
IF  (ZS.GT.Z)  Z=ZS 
CONTINUE 


S=0 

DO  3  K=1,M 
SS=0 

DO  4  J=1 ,  N 

4  SS=SS+ABS(  A(J,K) ) 
IF  (SS.GT.S)  S=SS 

3  CONTINUE 

Q=0 

DO  5  J=1 ,  N 
DO  5  K=1  ,M 

5  Q=Q+A( J,  K  )**2 
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Q=SQRT  (Q) 

AC=0. 

DO  6  J=1,N 

DO  6  K=1  ,M 

AB  =  ABS  ( A(  J,  K) ) 

IF  (AC.LT.AB)  AC=AB 

CONTINUE 

RETURN 

END 


